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ABSTRACT 

The Rossi X-ray Timing Explorer has detected nearly coherent oscillations in the tails of type I X-ray bursts 
from 17 low-mass X-ray binaries. The oscillations are thought to be generated by brightness fluctuations 
associated with a surface mode on the rotating neutron star The mechanism that drives the modes is, however, 
not understood, since the burning layer is stable to thermal perturbations. We show here via a linear perturbation 
analysis that, even under conditions when pure thermal perturbations are stable, nonradial surface modes may 
still be unstable by the e mechanism. Specifically, we find that, if helium-burning reactions supply a reasonable 
fraction of the outgoing flux during burst decay, nonradial surface modes will grow in time. On the other 
hand, the same modes are likely to be stable in the presence of hydrogen burning via the rp-process. The 
results naturally explain why oscillations in the decay phase of type I X-ray bursts are detected only from 
short-duration bursts. 

Subject headings: accretion, accretion disks — stars: neutron — X-rays; binaries — X-rays: bursts 



1. INTRODUCTION 



Type I X-ray bursts 
the surfaces of accreting 
stable hydrogen or helium 



are thermonuclear explosions on 
neutron stars, triggered by un- 
burn ing (jBabushkina et al.l 



[1971 iGrindlav & Heisa [TqtI iGrindlav et al.l 
Belian et'a l.1 Il976t l\yooslev & TaamI 119761: iJossI 
Maraschi & Cavalierd 119771: iLamb & LambI I1977L [ 
They typically have fast rise times of 1 s, exponential-like 
decays lasting ~ 10-100 s, and recurrence tim es of a few 
hours to days (for recent reviews, see ,Cummin3 l2004t 
IStrohmaver & Bildstenlllopeh . 

IStrohmayer et al.l (1 19961) were the first to detect coher- 
ent oscillations during a type I X-ray burst. Since their 
initial discovery, oscillations have been seen in lightcurves 
from helium-triggered bursts in 17 sou rces (Strohmayer & 
Bilds te n 2006 and refere n ces therein; iBhattacharvva et aU 
l2006t iKaaret et al.1 IIoOTL iBhattacharvval 120071) . as well 
as from one carbon-triggered s uperburst in 4U 1636-536 
(IStrohmayer & Markward t' 2002). Oscillations occur during 
the burst rise as well as the burst tail. The oscillation fre- 
quencies often increase by a few percent during the course 
of a burst, asymptoting to a value that is foun d to be sta- 
ble over years to within a few parts in lO-' (e.g.. iGiles et al.l 
I2OO2I: iMuno et aTl |2002|) . The asymptotic frequency is 
thought to correspond to the neu tron star spin frequency (e.g., 
IStrohmayer & Markwardtiri999l) . As of this writing, burst os- 
cillations have been detected only from sources with large 
values of the inferred accretion rate ((Muno et aLll2000l 120041 : 
Ivan Straaten et al.ll2001l:lFrancoll2001lK 

Oscillations during the burst rise are believed to be 
caused by the rotational mod ulation of a growing hot spot 
( IStrohmaver et al.lll997l 119981) . The thermonuclear instabil- 
ity that causes a burst is almos t certainly triggered initially at 
a point dJossI 1 978tlSharal 1 982h . and the resulting hot spot pre- 
sumably then expands and engulfs the entire s urface in ^ 1 s , 
the burst rise time (Frvxell & Wooslev 1982; iBildstenll 19951: 
IZingale et aT]|2001;,Spitkovskv et alji2002i) . At relatively low 
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accretion rates, ignition most likely occurs at the equator, and 
the resulting thermonuclear flame quickly spreads in longi- 
tude and generates an axisymmetric belt around the neutron 
star. Oscillations are not expec ted to be detected bec ause 
there is no azimuthal asymmetry dSpitkovskv et al.ir2002l) . At 
high accretion rates, however, ignition most likely occurs off 
the equator CCooper & Narayan .2007). The resulting ther- 
monuclear flame propa gates in longitude much more slowly 
dSpitkovskv et aljr2002l) . creating a relatively long-lived non- 
axisymmetric hot spot and leading to observable oscillations. 

While this simple model explains the oscillations observed 
during burst rise, it does not work for oscillations during 
burst decay since the data suggest that the entire neutron 
star surface radiates during this time. The most promis- 
ing explanation for oscillations in the burst tail is that they 
are generated by exci ted su r face modes on the neutron star 
dMcD ermott & Taam 19871 lEe^ 12001 iHeyll j200l l200l 
iLee & S trohmaver 2005; Piro & Bildsten 20(36^. In an im- 
portant paper. Hey], (.2004 ) identified surface r-modes as the 
most promising candidate. These modes travel backward (east 
to west) in the corotating frame and so the observed oscilla- 
tion frequency is reduced relative to the neutron star spin fre- 
quency. As the neutron star surface cools during the burst 
decay, the mode frequency decreases, thus explaining why 
the observed oscillation frequ ency increases. However, the 
expected frequency drift that Hevll d2004l) calculated is sig- 
nificantly larger than the observed drifts. |Piro & Bildst^ 
(|2p05b) showed that a surface mode i n the burning layer tran- 
sitions into a crustal interface mode dPiro & Bildstenll2005al) 
during the burst decay. This reduces the expected frequency 
drift to values that are in accord with the observed drifts. 

There is, however, one major unsolved problem: It is not 
clear what excites the surface mode during the burst decay 
and what keeps it active several seconds after the burst was 
initially triggered. It is highly unlikely that oscillations in the 
burst tail are produced by remnant asymmetry left over from 
the initial aspherical ignition. The cooling time at the burst 
peak is much less than the burst duration, so any initial asym- 
metry would have dissipated by this time. Also, since the 
modes identified by Heyl occupy only a small area of the neu- 
tron star surface around the equator, one requires huge oscilla- 
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tion amplitudes in order to explain the ^ 10% flux variations 
observed. It is hard to see how such a large amplitude could 
be maintained late in the burst if the modes were excited only 
at the time of ignition. 

In our view, a more reasonable explanation of the observa- 
tions is that the mode responsible for oscillations during burst 
decay is unstable. An instability would naturally explain the 
large oscillation amplitude. It would also explain why some 
type I X-ray bursts exhibit oscillations late in the burst decay, 
well after the time of the burst peak, when the thermonuclear 
flame has engulfed the entire stellar surface. 

In this investigation, we analyze the surface modes of a 
bursting neutron star, paying particular attention to the heat- 
ing and cooling processes in the burning layer We identify the 
conditions under which modes will grow and show that these 
conditions are likely to be satisfied for at least some type I 
X-ray bursts. Our results naturally explain why oscillations 
in the tails of bursts from nonmagnetic accreting neutron stars 
are excited specifically in short bursts. 

We begin in ^by deriving the equations that govern the 
mode dynamics. In ^we discuss the thermal structure of 
the accreted layer during the burst peak using a simple one- 
zone model. We combine the mode equations and the ther- 
mal structure model in ^ and analyze the resulting modes. 
We then determine the stability of these modes in ^ We 
derive the criterion for instability and demonstrate that sur- 
face modes are likely to be excited during the decay phase 
of type I X-ray bursts. We conclude in ^ with a discus- 
sion of the results. It should be noted that the mechanism 
for mode instability described here is different from the shear 
instability of zonal flows discussed by Gumming (20 p5j. I t 
is similar in spirit to the model of iPiro & BildstenI (120041) . 
who analyzed the stability of nonradial oscillations during 
thermally stable nuclear burning on the surface of a neutron 
star accreting at a super-Eddington rate, and to the models 
of McDermott & Taam ( 1987) and Strohmayer & Lee ( 1996), 
who found that the e mechanism could drive surface modes 
during thermally unstable nuclear burning. 

2. HEIGHT-INTEGRATED APPROXIMATION AND MODE 
DYNAMICS 

We make two principal simplifications in our analysis of 
the oscillation modes of the burning surface layer of a rapidly 
rotating accreting neutron star: 

(I) We assume that the mode is concentrated within a narrow 
equatorial belt. This allows us to ignore the spherical geome- 
try of the stellar surface and thereby obtain analytic solutions. 
As shown in the analysis of Longuet-Higgins (1968, see also 
Bildsten, Ushomirsky & Cutler 1996), the modes are indeed 
equatorially concentrated whenever the angular frequency of 
the star fl is much larger than the characteristic mode fre- 
quency. This condition is satisfied for rapidly rotating neutron 
stars with ~ few x lO"' radians per second. 

(II) We assume that each column of matter is always in radial 
hydrostatic equilibrium. This is a safe assumption since the 
vertical sound crossing time of the layer (^ 1 /is) is much 
shorter than the typical mode period (^ tens of ms). 

Making use of assumption I, we treat the equatorial belt of 
the star as a cylindrical surface of radius R and surface gravity 
g. We use coordinates x, y, z to represent linear displacements 
in (i) the azimuthal direction (east-west), (ii) the meridional 
direction (north-south) with y = Q located at the equator, and 
(iii) the radial direction (vertical) with z = located at the bot- 
tom of the burning layer, respectively. The coordinate x is 



periodic with a period of InR; by assumption I, the range of 
y of interest to us is small compared to R; since the layer is 
generally very thin 100 cm), the range of z is even smaller. 
Thus, over the volume of interest, the parameters R and g may 
be assumed constant. However, the Coriolis force, which de- 
pends on the projection of the stellar angular velocity vector 
on the local normal to the surface, does vary with position. 
Making use of assumption I, we write the Coriolis parameter 
/as 

(1) 



f = 2n-, 



which is valid for small excursions around the equator This 
is the well-known /3-plane approximation used in the study of 
Rossby waves in the atmosphere (e.g., Houghton 2002). 

In equilibrium, the density, temperature, and pressure of the 
layer are independent of x, y, and time t and are described by 
the functions po(z), Tq{z), po(z), respectively. Let us repre- 
sent the first-order perturbations of the density, temperature 
and pressure by p', T , and p', each of which is a function of 
X, y, z, and t. The equilibrium configuration has no motions 
as viewed in the rotating frame. The perturbations do involve 
motions, and we write the three components of the velocity 
perturbation as u, v, w. The continuity equation and the three 
components of the momentum equation then take the form 
dp' f du dv dw\ ^ 



dt 



dv dp' 

dw dp' , ^ 



(3) 



(4) 



(5) 



We now make use of assumption II and ignore the dynamics 
in the z-direction. That is, we assume that the pressure p at 
any given depth z in a local patch of the plasma is related 
simply to the local column density S above this z by 

p(x,y,z,t) = g / p(x,y,z,t)dz, (6) 

Jz 

where z = h corresponds to the surface of the layer. Let us 
define Eq to be the equilibrium column density of the burn- 
ing layer (a constant), E (x,y,t) to be the perturbed column 
density, and E to be the sum of these two. Integrating over z. 



E = / pdz = Eo + E 



(7) 



h = h^^ + h', (8) 
where, as before, the subscript corresponds to the equilib- 
rium value of a quantity and a prime corresponds to the linear 
perturbation. Similarly, we define a height-integrated pressure 



P = 



pdz = Po + P 



(9) 



We now integrate equations (l2]i-(|5]l over z to obtain the fol- 
lowing height-integrated dynamical equations: 



dt ^\dx dy 



= 0, 



„ du dP y 

^0^ + ^ 2r!Eo^v = 0, 

dt dx R 

dv dP V 

dt dy R 



(10) 



(11) 



(12) 
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where we have used equation ([T]i for /. In carrying out the 
height-integration, we have assumed that u and v are inde- 
pendent of z- This is a reasonable assumption so long as the 
mode under consideration has no nodes in the vertical direc- 
tion. Henceforth, our work is restricted exclusively to such 
modes (which are likely to be the most interesting ones for 
burst oscillations). 

We decompose the linear perturbations into modes in the 
usual way and write 



S {x,y,t) = E (y)exp(-iLUt + imx/R), 



(13) 



with similar expressions for P , T , u, v, etc. Here, uj is the 
mode frequency (as seen in the rotating frame), which is in 
general complex, and m is the azimuthal wave number, an 
integer. Substituting in equations (fT0t- (fT2l i. we obtain the 
following mode equations: 



iuj I im dv 

E + — u+ — = 0, 

Eo R dy 

im I in 

-iuju-\ P vv = 0, 

«Eo R 

1 dP in 

-icjv+ — — + —yu = 0. 
E() dy R 



(14) 
(15) 
(16) 



To close this system of equations we need an equation of 
state relating the effective two-dimensional height-integrated 
pressure perturbation P and the column density perturbation 
E . If we assume a simple adiabatic or polytropic equation 
of state, all the mode frequencies will be real and there will 
be no instability. Instability is possible only as a result of an 
interaction between the mode dynamics and thermal effects 
associated with nuclear reactions in the layer This motivates 
the discussion in the next section. 

3. THERMAL STRUCTURE OF THE BURNING LAYER 

For simplicity, we make use of a one-zone approxima- 
tion for the vertical structure of t he accreted layer (e.g. 
iFuiimotoet all [19811 iBildstenI [19981) . Thus, at each (x,y), 
we assume that the matter has constant density and write the 
height of the layer simply as 

ho=^, ^=-, h'^h-ho=^(^-^]. (17) 
Po p Pq yEo Po J 

We assume an ideal gas equation of state with adiabatic index 

k 



pkT T 
linip \p^ 



cv ■■ 



(j-l)pmp' 



(18) 



where s is the entropy per unit mass and cy is the specific 
heat per unit mass at constant volume. Since p = gE at the 
bottom of the layer (by the assumption of vertical hydrostatic 
equilibrium), we obtain by perturbing the equation of state 
that 

^4-^. (19) 
To Eo po 

Combining the one-zone approximation with the condition 
for hydrostatic equilibrium, we find for the layer height 



kT 
gpnip 



(20) 



and for the height-integrated pressure 



P= / pdz= / gp(h-z)dz=-gJ:h. (21) 
Jo Jo ^ 

Differentiating the latter, we obtain 

2 " ° Eo ho ) 2po Eo po ) ^ ^ 

To convert this to an effective equation of state of the form 
P = P (E ), we need a relation between p' and E . The rest 
of the section is devoted to deriving this relation. 

The entropy per unit volume of the matter in the layer 
evolves according to 



ds dF 



(23) 



where e is the energy generation rate per unit mass through 
nuclear reactions, and F is the energy flux. For the latter we 
use the radiative diffusion equation 

ac dT'^ 

F=^^^ (24) 
3k aE 

where k is the opacity of the plasma. Using the one-zone 
approximation and integrating over z, these relations become 

acT'^ 



■■^e-(F-Fb), 



(25) 



3kS ■ 

Here, Fb is the flux flowing into the base of the layer from the 
interior of the star, which we consider to be a constant, and F 
is the flux escaping from the upper surface of the layer 

We now consider linear perturbations of the entropy equa- 
tion. Let us begin with the entropy term. For small perturba- 
tions 

V_(7-iy' 

To Po 



■ Cv 



(26) 



Substituting for T /Tq from equation fT% and setting d/dt = 
-iuj, the perturbation of the entropy term in equation dZSl ) be- 
comes 



dt 



-iughoY^o 



IP 



(7-l)Eo (7-l)po 



(27) 



The nuclear energy generation rate e is in general a func- 
tion of both density and temperature. For small perturbations 
around the equilibrium, let us write this dependence via two 
indices 77 and i^. 



■eol — 
\Po 



where for later convenience we have defined eo = gho^h 
such that il^' is the characteristic heating time of the mat- 
ter through nuclear reactions. The perturbation of the heating 
term in the entropy equation dZSl ) is then 



(Ee) =T,oghonh 



(1/+I)— + (?7- 
^0 



^)^ 



Po 



(29) 



Similarly, for the radiative cooling term F, let us write the 
dependence of the opacity k on density and temperature as 

K{p,T) = Ko[ ^ \ (l^ , (30) 



KO — 

\Po 
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and also define fie in terms of the equilibrium escaping flux 

Fo, 

Fo = Sog/ioa-, (31) 

so that 51"' is the characteristic cooling time of the layer The 
perturbation of the cooling term then gives 



Equations (|39] l and ( l40t can be solved for E and u in terms 
of V and dv/dy: 



I 



(3-C)^-(4 + eC)- 
Eo po 



m^AS) 
i 



m y 



IVLR R 



A dv 

md^ 



(32) 



(X^-m^AS) 



\-v-mA5R — 



R 



dy 



(42) 



(43) 



Substituting equations (|27] i. ( |29l ) and ( [32] i in the perturbed 
entropy equation 



Substituting these in equation (|40] |. we obtain the following 
second-order differential equation for v: 



= (Se) -F 



(33) 



dy^' 



A 



= 0. 



(44) 



and rearranging terms, we obtain an expression for p'/A) 
terms of E /Eq. Substituting this in equation ( l22l i we finally 
obtain the effective equation of state of the matter in the layer. 
The result is 

p'=Agho^', (34) 
where A is the following dimensionless quantity, 

27 

1 + j[(7- l)/(27- 1)] [(^^A/^)(2r?- ^+ 1) + (0,/c^)(5 + 2^- 0] 



1 + i[(7- l)/7] [(^^a/^)('7- '^) + {n,/Lo){4 + ^-0] 



(35) 



For a given mode, A is a constant; it is, in general, complex. 

For pure adiabatic perturbations (i.e., no heating or cool- 
ing), A simplifies to (27-l)/27 and the height-integrated 
pressure perturbation satisfies 



Equation ( l44b may be directly compared to equation (8.2) in 
Longuet-Higgins ( 1968), noting that our m and AS are called 
s and 1 /e in tha t paper. The term -m^ in our equation is not 
present in Long uet-Higginsl (ll968h: however, this term is neg- 
ligible whenever the assumption of an equatorially concen- 
trated mode (on which our analysis is based) is valid. 

Before proceeding to solve the mode equations, we note 
that we are not interested in modes with m = 0, since they will 
not produce any observable oscillations. Nor are we interested 
in modes with very large values of \m\ or modes with many 
nodes in the y-direction, since these will again have little de- 
tectable variations in the observed flux. In the following, we 
limit our discussion to the lowest order eigenfunctions in the 
y-direction. The generalization to higher order modes is, of 
course, trivial. 

Equation (|44] | is of standard form and its solutions are Her- 
mite functions. The lowest order eigenfunction is given by 



v(y) = ICVlRS^I^ exp 



P 



(27-l)S 



y 



2(A(5) 1/27^2 



(45) 



(36) 



7 So 

Thus, the two-dimensional fluid behaves as if it has an effec- 
tive adiabatic index 

(27-1) 



7eff = 



7 



(37) 



where C is an arbitrary dimensionless constant that mea- 
sures the mode amplitude, and we have introduced the factor 
iriRS^I^ to scale the velocity by the sound speed. The func- 
tion ( |45l l is a solution of equation ( l44l l provided the eigenvalue 
A satisfies the following quadratic equation 



If the fluid is incompressible in three dimensions (7 00), for 
example, the height-integrated two-dimensional system has 
an effective index 7eff = 2. By noting this correspondence, 
one can easily map many of the equations in the present paper 
to corresponding results in Longuet-Higgins 's (1968) analysis 
of incompressible ocean waves. 

4. ANALYSIS OF MODES 

Let us define the following two dimensionless quantities: 



\^-m{A6)^l^\-(A6fl^ = Q. 



(46) 



A = 



5 = 



(38) 



2r2' " A9?R'^' 
Substituting equation ( |34l ) in equations (fl4]l-(fT6]l and rewrit- 
ing in terms of A and 5, we obtain the following set of equa 
tions for E , u and v: 



E im 1 dv 
~'^%^7nR"^7hTy^^' 

E' y 

-/Am + 2imQ.RA6— - - v = 0, 
Eo R 



-iXv + inR^AS- 



dy \ Eo 



+ -u = 0. 



(39) 
(40) 
(41) 



This gives two roots, which correspond to the g-modes dis- 
cussed in § 4.1. Actually, there is a third root, A = -miASy^-^, 
but it is spurious since it causes the denominators in equa- 
tions (l42l i and ( |43T l to vanish. This spurious root is replaced 
by another root, which we discuss in § 4.3. 

In order for the solution ( |45] | to be physically valid it must 
decay rapidly with increasing \y\. This means that, in taking 
the square root (A(5)'/2, we must select the solution with a pos- 
itive real part. Further, in deriving the equations, we assumed 
that the mode is concentrated around the equator. For this to 
be true, we see from equation ( l45T l that we require the modu- 
lus of (A(5)'''2 to be <C 1. The quantity A 1, so we require 
5'/2 <^ 1 poi- a typical neutron star (M 1.4M0, 10 
km), ^ - 2 X 10'*cms-2, ho - 200 cm and Q/ln - 500 Hz. 
We then find that S^^^ ^ 10"' which is indeed small. The 
mode is thus restricted to \y\/R < 0.2, which is sufficiently 
concentrated around the equator for our approximations to be 
valid. 

The next-order solution of equation (l44l i is 



v(y) = 2CnRS^^^ 



{ASy/'^R 



exp 



2(A5)'/2/;2 



(47) 
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We have included an extra factor of (AJ)"''''* to allow for the 
typical value of \y\/R, and so C is again a measure of the di- 
mensionless mode amplitude. This solution gives the follow- 
ing condition on the eigenvalue A: 



as' 



m 

a" 



:0. 



(48) 



Longuet-Higgins (1968, see also Heyl 2004) has shown 
that the surface modes of a rotating sphere separate into three 
kinds: ^-modes, r-modes and Kelvin modes. We discuss each 
of these below, focusing on the lowest order eigenfunction in 
each case. 

4.1. g-Modes 

The lowest order ^-modes correspond to the solution ( l45l l. 
with the dimensionless mode frequency A obtained by solving 
equation (|46] |, 

A « ±{A5yi\ (49) 

where we have made use of the fact that 5 is small. For a given 
m there are two modes, one of which moves towards the east 
(the mode with the + sign) and the other moves towards the 
west (the mode with the - sign). Using the typical numerical 
values given above and setting A 1, we find that the typi- 
cal mode frequency is ^ 200 Hz, which is fairly large. The 
eigenfunction has the form (we show only the leading terms) 

„2 



S(y) 

So 
u(y) 



y 



'.±iC 



iC 

Ai72 {ASyl^R 

y 



exp 



y 



(ASyl^'R 



exp 



2(A5) 1/2/^2 



2(A5)i/27;2 



: Cexp 



y^ 



2(A5)i/2/;2 



(50) 



(51) 



(52) 



The g-modes have a lower growth rate than the other two 
modes discussed below. For this reason, and also because 
they have larger frequencies, these modes are probably not of 
great interest for burst oscillations. 

4.2. r-Modes 

The lowest order r-modes are given by the solution (l47T i. 
with the eigenvalue obtained by solving the cubic dispersion 
relation ( |45] l. Two of the roots correspond to higher order g- 
modes and we ignore them. The third root gives the r-mode. 
In the limit {A5y/- <C 1, this root is approximately equal to 

X^-lmiASyl^. (53) 

For m= 1, 2, the typical mode frequencies are ~ 10, 20 Hz, 
respectively. The negative sign on A indicates that the r-mode 
always travels towards the west. Both of these properti es are 
Just w hat are required for explaining burst oscillations (iHeyll 
[2001 . 

The eigenfunction of the r-mode is given by 



S(y) 



u{y) 
v{y) 



3iC 



8mA^/HA5y/* 



1 + 



2/ 



(ASy/^R^ 



exp 



y 



4.3. Kelvin Modes 

The Kelvin mode has v identically equal to zero, and so it 
cannot be obtained by solving the differential equation (l44l i. 
Instead, one has to go back to the primitive equations ([39l)- 
(l40l i. The eigenvalue of this mode is 

A = m(A5)'/^ (57) 
and the two nonzero components of the eigenfunction are 

^'(y) c 



So 



exp 



= Cexp 



2(A<5)V2/;2 

y' 



liASy^R^ 



(58) 



(59) 



The Kelvin mode always moves towards the east. For a 
given m, it has a frequency equal to three times that of the r- 
mode. The motions of the plasma are entirely in the east-west 
direction (v = 0) whereas they occur in both the east-west and 
north-south directions in the r-mode. 

5. GROWTH RATE OF MODES 

Equations (l49T l, ( l53T l and (l57T i give the frequencies of the 
lowest order g-, r- and Kelvin modes. If A is real, then all the 
frequencies are real and there is no mode growth. However, 
we showed in § 3 that A is in general complex. Therefore, the 
modes will either grow or decay with time. 

5.1. Criterion for Mode Growth 

Let us focus on the r-mode and Kelvin mode for now. Re- 
call that these modes have angular frequencies uo ^ \Q- radi- 
ans per second. In contrast, the cooling time of the accreted 
layer is of order a second, so that the parameter fif ~ 1 s"' . 
The parameter 57^ is even smaller during the decline of a burst. 
We are thus permitted to treat Q,c/lo and Vlk/to in equation 
(|35] | as small quantities. We can then Taylor expand ( |35] | to 
obtain 



1 1/2 



,1/2 



1 + / 



-a 



(60) 



where 



,1/2. 



(27-1) 
27 

(7-1) 
27(27- 1) 

(7-1) 
27(27- 1) 



1/2 



(61) 

[7 + r? + (7-l)i^], (62) 

[37-4-C-(7-l)C]. (63) 

1/2 



As discussed earlier, we should pick the positive root for A^' 
to ensure that the eigenfunction decays at large y. 

Substituting (|60] | in equation i53[ . we find for the frequency 
of the r-mode 



2mVL 



{Ar6) 



1/2 



l+i a p — 



(64) 



2(A5)'/2/?2 

(54) 



9iC 



= C 



m(A^)i/4 

y 



2/ 



2QR6^/^ (ASy/'^R 



exp 



3(A(5)i/2/?2 

,2 



exp 



r 



2(A(5) 1/27^2 



,(55) 



Since the magnitude of the imaginary term on the right is very 
small, we immediately obtain 

^■^{AKSy/^ + i{an,-pn,). (65) 



y 



2(A5)i/2/?2 



(56) 



In a similar manner, the frequency of the Kelvin mode is 
found by substituting ( l60b in ( 157b : 



Lo « 2mn{AR5yl^ + i{anh-l3nc) 



(66) 
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In the case of the ^-mode, ( |49] l shows that A oc {ASy^"* rather 
than (A(5)'/^. Consequently, the imaginary part of uj is reduced 
by a factor of 2: 

UJ « ±2n(ARSy^^+^(aQh-f3nc). (67) 

A mode is unstable if its uj has a positive imaginary part. 
Thus, all three modes give the same condition for instability, 
viz., 

anh-f3n,>Q, (68) 

with a and (3 defined in equations ( |62] ) and ( 1631 ). Before 
applying this result to burst oscillations, we first discuss the 
physics of the instability. 

5.2. Physics of the Instability 
5.2.1. Local Thermal Stability 

We begin by deriving the standard criterion for thermal sta- 
bility of a burning layer We consider linear perturbations of 
the entropy equation dZSl ). but we assume that there are no 
nonradial perturbations, i.e., S = 0. From equations (fT9] l. 
dlTli, (lig, and setting E' = 0, we find 



d\nT 



(7-1) dt 



(69) 



Clearly, the accreted layer is thermally unstable when the right 
side of the equation is positive. Thus, the criterion for thermal 
instabihty is 

§i>ii«, (70, 

\lc v — rj 

Consider the accreted layer on a neutron star prior to the 
onset of a type I X-ray burst. The layer undergoes steady 
thermonuclear burning at some rate. Steady state requires in 
general Soeo = pQ-Ft. However, the flux Pt, entering the base 
of the accreted layer is typical ly much less than th e flux Pq 
escaping from the surface (e.g.. lBrownll2000i |2004|) . so Pq ^ 
Pb- Therefore, we expect Soco ~ Pq, which is equivalent to 
the condition 51/, « flc- The criterion for thermal instability 
then becomes 

j.-^>4 + e-C, (71) 



which is equivalent to equation (23) of BildstenI (119981) . 

The thermal instability is driven by variations in T. Con- 
sider a positive temperature perturbation at constant pressure 
(the pressure is constant because we have assumed hydrostatic 
equilibrium and a constant S). This leads to a corresponding 
decrease in p according to equation (fT9l l. If the marginal in- 
crease in the nuclear energy generation rate e (xi v-rf) ex- 
ceeds the marginal increase in the effective cooling rate of the 
layer P /E oc (4 + ^- Q, the temperature of the accreted layer 
increases further and a thermonuclear runaway ensues. This 
explains the form of the two factors on either side of equation 
( TTTT i as well as the direction of the inequality. 

The thermal instability is paramount for understanding 
the triggering of type I X-ray burs ts, and there is consid- 
erable literature on this subject ('Schwarzschild & Harm] 
1965t [Hansen & van Horiil i l975: Fuiimoto et al. 19811 
Paczvnskil 11983a'; 'Fushiki & Lamb 1987; Bildsten 1991 
Narayan &Heyl 2003; Cooper & Narayan 2006^3 
Fisker et al.l 120061 |2007). However, once a burst has 
been initiated, the temperature of the burning matter in- 
creases dramatically until the burning layer finds a new 



quasi-steady-state in which the temperature is typically 
~ 10^ K. In this new state, the indices v, rj, ^ and ( have 
different values. In particular, the temperature index i' of the 
nuclear energy generation rate e, which is large at the cooler 
temperatures characteristic of burst ignition, now becomes 
much more modest. As a result, the condition ( ItTI i is no 
longer satisfied, and the layer is thermally stable. Thus, as 
far as pure thermal perturbations are concerned, the burning 
layer during a burst is stable and the observed oscillations in 
the burst tail cannot be attributed to thermal instability. 

5.2.2. Nonradial Mode Stability 

The instability that we have analyzed in this paper is as- 
sociated with nonradial surface modes, and it depends cru- 
cially on the fact that S is perturbed. The physics is similar 
to that which causes no nradial stel lar pulsations in Cepheids 
and other variable stars (ICoxlll980l) . The one big difference is 
that the primary driving mechanism here is nuclear reactions 
(the so-c alled e mechanism) r a ther th an opacity ( the k mech- 
anism). 'McDe rmott & TaamI (Il987h and IStrohmayer & Led 
tl996) investigated the driving of surface modes via the e 
mechanism in thermally unstable stellar envelopes. In con- 
trast, we have investigated the driving of surfa ce modes in 
thermally stable envelopes. As iPiro & BildstenI (|2004) have 
emphasized, the stability criterion for nonradial perturbations 
is distinct from the criterion for thermal stability. 

In stellar pulsation theory, the following simple criterion is 
used to decide whether a particular layer of matter is a driving 
or damping region for pulsations (Cox 1980): A region that 
gains heat when compressed is a driving region, whereas a 
region that loses heat when compressed is a damping region. 
In the nonradial surface modes of our problem, compression 
is induced by variations in E. Since we are considering a 
simple one-zone vertical structure for the layer, we merely 
have to determine whether the specific entropy of the plasma 
increases or decreases when S increases. If the entropy in- 
creases with increasing E it means that the plasma gains heat 
when it is compressed and the mode is unstable; if the entropy 
decreases, the mode is stable. 

In equilibrium, an unperturbed column satisfies Eoeo = 
(Po-Pb)- Let us imagine compressing the column laterally 
so that the surface density increases to Eq + S . Let the com- 
pressed column relax to a new thermal equilibrium state. As- 
suming Pt, is unchanged, we require (Ee) = P' . Substituting 
for these quantities from equations ( l29b and ( |32] |. we can solve 
for p' / po in terms of E /Eq: 



P_ 

Po 



-(t/+l)f7,, + (3-C)flc E 



(72) 



(7j-iy)nh + i4 + ^-Onc Eo 
From equation (T% . we can also obtain an expression for 
r'/Ti) in terms ofE'/E,,: 



T' _ (ry+l)»„ + (e+l)», E 

7b (Tj-u)nh+(4+^-on,^o' 



(73) 



We may now use equation (1261 ) to estimate the perturbation in 
the specific entropy of the plasma: 



[7 + 77+(7-l)i/]r!A-[37-4-e-(7-l)C]f^c S 

Cv 



i7j-iy)n„ + i4 + ^-On, 



So' 



(74) 



By the simple criterion mentioned earlier, the layer of 
plasma under consideration will be a driving layer for pulsa- 
tions if s' increases with increasing E , i.e., if the coefficient 
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on the right side of equation (|74] | is positive. During the burst, 
the layer is thermally stable, and so by equation (fTOl l the de- 
nominator of this coefficient is positive. Thus, for nonradial 
instability, we require the numerator also to be positive, i.e., 

[7 + 7y+(7-l)i/]a-[37-4-e-(7-l)C]a.>0, (75) 

which is identical to equation ( l68T l. 

The analysis presented here assumes that the column comes 
into thermal equilibrium after it is compressed. This is not a 
good approximation in the case of a mode, especially if the 
mode frequency is much larger than the heating/cooling rate, 
as in our problem. Nevertheless, it is easy to see that the same 
criterion for instability must hold. Consider a column of mat- 
ter that undergoes a positive S perturbation as a result of lat- 
eral compression induced by an oscillating mode. Since the 
compressed column is not in thermal equilibrium, its specific 
entropy will begin to increase. Of course, the entropy will not 
have time to reach its equilibrium value since the mode varies 
rapidly. Nevertheless, by the time the compression begins to 
reduce, the entropy of the plasma will be a little larger than 
before the column was compressed. The increased entropy 
will be reflected in an increased effective pressure P, and so 
the expansion will be a little more energetic compared to the 
initial compression. As a result, the amplitude of the next 
half-cycle of the wave will be a little larger The process is 
obviously self-sustaining and the wave amplitude will grow 
exponentially with time. 

5.3. Application to Type I X-ray Bursts 

Surface modes during a type I X-ray burst will be unstable 
whenever the condition ( |68] | or ( iTSl l is satisfied. At the high 
temperatures that are present during a burst, the op acity scales 
according to the approximate formulae given in iPaczvriskil 
( Il983bl) . which give -0.1 < C < 0, -0.5 < C < 0. We assume 
^ = 0, C = -0.25. During the decline phase of a burst, the es- 
caping flux is obviously larger than the rate at which heat is 
added via nuclear reactions, i.e., il^ > fi/,. The ratio of the 
two rates, however, depends on details which are beyond the 
scope of this work. 

When gas pressure dominates, we have 7 = 5/3, and the 
condition for nonradial instability becomes 

^ > ^ r76) 

10 + 6?7+4z^' 

On the other hand, when radiation pressure dominates, 7 = 
4/3, and the condition is 

^> ^ . (77) 

Vic l6+l2ri+4iy 

Whether or not these conditions are satisfied depends on the 
kind of burning that takes place during the burst. 

5.3.1. Helium Burning 

The nuclear flow during a helium-rich burst includes few 
weak interactions, and thus nuclear burning is rapid and oc- 
curs over a relatively short timescale. Consequently, Vti, is 
likely to be a considerable fraction of ftc during the burst. 
Furthermore, if heating during the burst is primarily through 
helium-burning reactions, we expect e to depend on both den- 
sity and temperature. Depending on whether the triple-a re- 
action or some other (two-particle) a chain reaction is impor- 
tant, ry = 2 or 1 . As for i^, we expect it to be of order a few. 

The temperatures achieved during a burst can often exceed 
10^ K, in which case radiation pressure becomes important 



( lJosslll977l: IWooslev et"ani2004 . This is especiafly true of 
really powerful, e.g., photospheric radius expansion (PRE), 
bursts. If radiation pressure dominates, the instability crite- 
rion is given by ( iTTl i. which is easily satisfied even if nuclear 
burning supplies only a few percent of the cooling flux. Thus, 
instability is very likely in such bursts. Observationally, a 
good number of bursts with osciUations are indeed fo und to 
be PRE bursts dMuno et al.ll2000HGallowav et al.ll2007l) . PRE 
bursts are typically helium-rich, and thus we conclude that 
helium burning is conducive to burst oscillations. 

If gas pressure dominates, then the instability criterion is 
given by ( |76] l, which for the typical values of 77 and ly given 
above, will be satisfied so long as f2/, is greater than or of order 
ric/4. The exact amount of nuclear burning during the decay 
phase of a helium burst is unclear. This issue is complicated 
by the fact that, although observations indicate that short, 
helium-dominated bursts occur at high accretion rates, all cur- 
rent detailed numerical burst models predict t hat such bursts 
occur only at low accretion rates. However. IWooslev et al.l 
(120041) show that nuclear burning is substantial for the major- 
ity of the duration of these bursts. Thus, it is probably not 
unreasonable to expect helium-rich bursts to have oscillations 
even when they are gas-pressure dominated. 

Helium burning is typically rapid, and so bursts in which 
helium burning dominates are relatively short, no more than 
several seconds. Hence we suggest that short bright bursts are 
the best places to look for burst oscillations. 

Given the short duration of these bursts, one issue is 
whether there will be enough time for the mode to grow to 
a nonlinear amplitude. The growth rate is given by equa- 
tion (|67] | and is generally a fraction of the heating/cooling rate 
(since a, f3 are typically less than unity). The burst duration, 
on the other hand, is a factor of several longer than the heat- 
ing/cooling time scale. Thus the growth time is expected to be 
comparable to the burst duration. It is then a delicate question 
whether or not the oscillations will have enough time to grow 
to a large amplitude before the burst dies away. This can be 
decided only with detailed calculations. 

5.3.2. Hydrogen Burning 

In hydrogen-rich bursts, much of the h eating is supplied by 
hydrog en burning via the rp-process of [Wallace & WooslevI 
(119811) . The nuclear energy generation rate is limited by a 
long series of slow /3-decays which act as bottlenecks in the 
nuclear flow. Because of this, the reaction rate per unit mass 
is less sensitive to density and temperature, i.e., rj and i' are 
closer to zero, compared to the case of helium burning. This 
means that flh / ilc has to be larger than in the case of helium- 
burning bursts before the instability will be triggered. Fur- 
thermore, because of the slow r/?-process, these bursts are not 
very bright at the peak and are unlikely to become very radi- 
ation dominated. As a result, they need to satisfy the more 
stringent criterion ( |76] | for instability. For all these reasons, 
we do not expect these bursts to exhibit oscillations. 

Hydrogen-rich bursts usually last substantially longer than 
helium-burning bursts because of the slow r/?-process. Ob- 
servationally, we thus suggest that long-duration bursts will 
in general show fewer episodes of oscillations. 

5.3.3. Oscillations During Burst Rise 

As explained in §1, oscillations during the burst rise are 
interpreted in terms of a non-axisymmetric burning spot asso- 
ciated with the initial burst trigger. However, unstable nonra- 
dial modes of the kind discussed in this paper should also be 
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considered. During the rise, cle arly fli, > il c, and so mode in- 
stability is quite likely. Indeed. iMuno et al] ((2002) found that 
bursts from which oscillations are detected during the rising 
phase show the largest frequency drifts, and they concluded 
that the mechanism that generates oscillations in the burst tail 
begins with the initial burst trigge r. Perhaps one could int er- 
pret the propagating spot model of lSpitkovskv et al] (l2002h as 
a nonlinear traveling mode. On the other hand, the burst rise 
is rather rapid, and it is possible that there is not enough time 
for a mode to grow to saturation amplitude. 

5.3.4. Amplitude of the Oscillations 

An unstable mode is likely to saturate at a dimensionless 
amplitude C < 1. Thus we expect E /So < 1 and, by equa- 
tion (|32] |. the amplitude of variation of the local escaping 
flux F to also be of this order. In the case of /? Cephei 
stars, temperature fluctuations in nonradial pulsational modes 
are typically about ±10% (Berdyugina et al. 2003), which 
correspond to flux variations of ±50%. Let us assume a 
similar amplitude for F in our problem. Since the nonra- 
dial modes we have considered occupy only about 20% of 
the surface of the star (see the discussion in ®, the ob- 
served flux variations in burst oscillations are expected to be 
^ 10%. This estimate is in accord with the observed ampli- 
tudes, which are < 15% (M uno. Ozel. & Chakrabarty„2002l: 
IStrohmayer & Bildstenll2006h . 

5.3.5. Radial Structure of the Mode 

The vertical thickness of the burning layer varies as in 
equation STJl . Thus h' /ho is of order S'/Eq, which itself is 
roughly the transverse displacement divided by the the wave- 
length of the mode. The latter quantity is of order h^/R, so 
we expect the ratio of the radial displacement to the trans- 
verse displacement to be of order Iiq/R ^ 10 "^. This value is 
consistent with the equivalent ratio f^/^j^ of iPiro & BildstenI 
(|2004). Figure 4 of lPiro & BildstenI (120041) illustrates that 
is roughly constant within the burning layer and drops rapidly 
below the burning layer In our simple model, ^± is in effect 
assumed to be constant throughout the burning layer and to 
go to zero immediately below the layer Our approximation 
is thus reasonab l y con sistent with the more detailed model of 
iPiro & Bildst"enl (|2004 . 

5.3.6. Frequency Drift 

A general feature of burst oscillations is that the frequency 
drifts upward during the course of a burst, typically by about 1 
to 4 Hz (e.g., Muno et al. 2003; Piro & Bildsten 2005b). Heyl 
(2004) argued that this is the result of the mode frequency in 
the rotating frame (the real part of uf) decreasing during the 
course of the burst as the plasma cools. For a backward mov- 
ing mode such as the r-mode, there would be a corresponding 
increase in the observed oscillation frequency. 

For an m = 1 r-mode, the typical value of Re(a;/27r) ~ 
I0(r/10''K)i/2 Hz (Piro & Bildsten 2005b). As a burst de- 
cays from its peak, the temperature is expected to vary by 
about a factor of 2, from say lO** K to 5 x 10^ K, and the cor- 
responding change in Re (w /2tt) is expected to be about 3 Hz. 
This estimate is comparable to the upper end of the range of 
observed frequency drifts, but is perhaps a little too large to 
be consistent with the smaller drifts 1 Hz seen in some sys- 
tems. To solve this problem (though we are not convinced that 
there is indeed a serious problem), Heyl (2004) suggested that 



the relevant modes may be some kind of photospheric distur- 
bances, while Piro & Bildsten (2005b) proposed a transition 
from a surface wave (of the sort described in the present pa- 
per) to a crustal interface mode. 

6. DISCUSSION 

As mentioned in ^ oscillations during the decay phase of 
type I X-ray bursts are detected predominantly when the in- 
ferred accretion rate is high. In ^5.31 we deduced that oscil- 
lations should occur only in short, helium-dominated bursts. 
Observations indicate that the burst decay timescale decreases 
with i ncreasing M dvan Paradiis et al.lll988t ICornelisse et al.l 
l2003t iGalloway et al.ll2007l) . Therefore, we suggest that os- 
cillations are detected predominantly at high accretion rates 
primarily because the burst durations are lower at these rates. 
To our knowledge, all bursts from neutron stars that are not 
accretion-powered pulsars and that exhibit oscillations in their 
tails have short durations. 

We stress that our work may not apply to accretion- 
powered millisecond pulsars from which burst oscillations 
have been detected. There are several differences between 
bursts from pulsars and those from non-pu lsars (see, e.g., 
iPiro & BildstenI I2005bt IWatts & Strohmayeri i2006). In par- 
ticular, oscillations have been detected only in long bursts 
from the two pulsa rs that ex hibit burst oscillations, XTE 
J18I4-338 ( StrohmaYereLaLl|2003) and SAX J1808.4-3658 
(lin'tZand et al.l l2001|) . Thus, we suspect that a different 
mechanism generates burst o scillations in accretio n-powered 
millisecond pulsars (see, e.g.. lLovelace et al]|2007l) . 

There are several issues that we have not addressed. 

First, we have not attempted to calculate the crucial heating- 
to-cooling ratio, which figures prominently in our 
instability criterion (|68] |. dTSl l. This can be done only with 
detailed multi-zone time-depend ent burst simulations (e.g., 
IWoosley et ani2004tlFisker et al.l l2006). 

Second, we have only carried out a linear stability analysis, 
and we had to use a plausibility argument to estimate the mode 
amplitude at saturation and the amplitude of the observed flux 
variations. To obtain a more accurate estimate, one needs to 
carry out full nonlin ear calculations, perhaps along the lines 
of Spit kovsky et"an (B002). 

Third, it appears from equations ( |65] | and (l66l l that the r- 
mode and Kelvin mode have identical growth rates and should 
therefore both be observed. Also, the linear growth rate does 
not depend on the azimuthal mode number m. However, the 
observ ed oscillati ons seem to consist predominantly of the r- 
mode ( Hevll2004l) with m=l. Differences between the modes 
might emerge in the nonlinear regime, driven perhaps by the 
fact that the frequency of the r-mode is a third that of the 
Kelvin mode and the spatial scale of the m = 1 mode is larger 
than that of higher m modes, but it is mere speculation at this 
point. 

Finally, oscillations have been detected in the decay phase 
of one superburst. We suspect that a mechanism similar to the 
one we have described for helium-burning type I X-ray bursts 
generates these oscillations as well. However, unlike normal 
helium-triggered bursts, the physics of the ignition and subse- 
quent nuclear burning of superbursts are not well understood 
(see iWeinberg et al.l20 06. for some recent work on this topic), 
so we cannot make a definitive statement regarding superburst 
oscillations. 

We would like to thank Tony Piro for his critical reading of 
an earlier draft of this manuscript, Lars Bildsten for helpful 
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